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Thermoelectric  devices  (TEs)  can  achieve  direct  conversion  of  heat  and  electricity  by  semiconductor 
materials,  coupling  of  heat  transfer  and  electric  conduction  is  important  to  accurately  predict  the 
performance  of  TEs.  This  paper  develops  a  general,  three-dimensional  numerical  model  of  TEs  with 
consideration  of  coupling  of  temperature  field  and  electric  potential  field.  The  model  is  used  to  figure  out 
the  performance  of  thermoelectric  coolers  (TECs)  with  the  temperature-dependent  thermal  conductivity, 
electric  conductivity,  and  Seebeck  coefficient  of  semiconductor  materials.  A  miniature  TEC  is  considered 
and  Bi2(Teo.94Se0.o6)3  and  (Bi0.25Sbo.75)Te3  are  selected  as  the  n-type  and  p-type  thermoelectric  materials, 
respectively.  The  effect  of  parameters  such  as  the  temperature  difference  and  the  current  is  investigated 
under  conditions  of  variable  material  properties  as  well  as  radiation  and  convection  heat  transfer 
occurred  between  the  TEC  and  the  ambient  gas.  The  results  show  that  the  variable  properties  and  the 
heat  losses  to  the  ambient  gas  have  significant  effects  on  the  cooling  capacity  and  the  coefficient  of 
performance  (COP)  of  the  TEC.  Three-dimensional  temperature  distributions  within  the  semiconductors 
is  observed  under  convective  boundary  condition  and  it  becomes  remarkable  at  large  temperature 
differences  and  high  currents. 

©  2012  Elsevier  Ltd.  All  rights  reserved. 


1.  Introduction 

Recently,  the  research  and  development  of  thermoelectric 
devices  (TEs)  have  attracted  a  great  deal  of  attention  because  of 
their  potential  applications  in  green  energy  and  energy  manage¬ 
ment  [1—3],  TEs  can  be  classified  into  two  different  groups;  one  is 
the  thermoelectric  generator  (TEG)  which  converts  heat  into  elec¬ 
tricity  [4—15]  and  the  other  is  the  thermoelectric  cooler  (TEC) 
which  converts  electricity  into  heat  [16—43].  Compared  to 
conventional  mechanical  providers  of  heating,  cooling  and  electric 
power  generation,  such  as  air  conditioners,  refrigerators,  heat 
pumps  and  turbine  engines,  TEs  are  more  robust,  compact  and 
noiseless  because  they  are  solid  state  heat  engines  that  use 
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electrons  and  holes  as  working  fluid,  so  that  neither  moving  parts 
nor  real  fluids  are  needed  [1,2], 

TEs  contain  over  hundred  pairs  of  thermoelectric  elements  (p-n 
junctions)  connected  electrically  in  series,  but  thermally  in  parallel 
between  two  planar  substrates.  The  TE  performance  is  closely 
related  to  the  figure-of-merit  of  thermoelectric  materials,  ZT,  which 
is  defined  as  ZT  =  cwrT/A,  where  a  is  the  Seebeck  coefficient,  a  is  the 
electric  conductivity,  X  is  the  thermal  conductivity,  and  T  is  the 
absolute  temperature  at  which  the  properties  are  measured  [2,3]. 
The  materials  with  large  ZT  can  improve  the  TE  performance 
significantly  [1—3],  thus,  high  electric  conductivity,  high  Seebeck 
coefficient,  and  low  thermal  conductivity  are  desired  for  thermo¬ 
electric  materials.  In  addition,  the  TE  performance  can  also  be 
improved  by  optimization  of  its  geometric  design  and  operating 
conditions. 

So  far,  a  number  of  models  have  been  proposed  to  analyze  the 
heat  transfer  and  performance  of  TEs.  The  models  can  be  divided 
into  two  types.  One  is  the  energy  equilibrium  model  or  zero¬ 
dimensional  model  [4-8,11,18,21,22,24,25,27,29,31,33,39],  where 
heat  balance  equations  (algebraic  equation)  at  the  hot  junction  and 
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the  cold  junction  are  constructed,  respectively,  without  solving 
energy  differential  equation.  The  advantage  of  these  models  is  that 
it  can  obtain  the  analytical  expressions  for  the  TE  performance,  but 
their  accuracy  is  limited  due  to  gross  simplifying  assumptions.  For 
example,  Joule  heat  and  Thomson  heat  are  generally  assumed  to  be 
equally  distributed  to  the  hot  and  cold  junctions.  The  other  is  the 
one-dimensional  [9,17,19,23,26,28]  or  three-dimensional 

[10.30.38.43]  heat  transfer  model,  where  Fourier’s  heat  conduc¬ 
tion  equation  (differential  equation)  with  Joule  heat  and/or 
Thomson  heat  as  inner  heat  sources  is  modeled  and  solved.  Then 
the  TE  performance  can  be  obtained  by  solving  temperature 
distributions.  In  these  two  types  of  models,  the  p-n  junction  is 
usually  simplified  as  an  individual  bulk  phase  so  that  the  difference 
in  the  thermal  behavior  between  p-type  and  n-type  semi¬ 
conductors  is  not  possible  to  be  distinguished  and  evaluated  [4— 
11,17—19,21—29,31,33,39].  Recently,  some  researches  treat  the  p- 
type  and  n-type  semiconductors  as  two  separate  parts  and  consider 
the  difference  in  temperature  between  the  two  elements 

[30.38.43] ,  Furthermore,  in  these  models,  the  current  is  assumed  to 
be  transported  only  along  the  height  direction  (z  direction)  of 
semiconductors  with  uniform  current  density  over  the  cross- 
section  (x — y  plane)  of  semiconductors,  thus,  the  current  density 
vector  J  ( x,y,z )  within  the  semiconductors  can  be  simplified  as 
a  scalar,  or ,JX  =Jy  =  0,Jz  =  I /A,  where  1  is  the  total  current  through 
the  TE  element,  and  A  is  the  cross-sectional  area  of  semiconductors. 
So,  the  electric  potential  equation  is  not  needed  to  be  solved. 

In  fact,  there  is  coupling  of  temperature  field  and  electric 
potential  field  within  semiconductors.  First,  the  electric  potential 
through  semiconductors  includes  two  parts:  the  electric  potential 
due  to  current  flowing  through  the  semiconductors  and  the  ther¬ 
moelectric  potential  produced  by  temperature  gradient  due  to 
Seebeck  effect.  The  Seebeck  coefficient  and  electric  conductivity  of 
semiconductor  materials  are  all  temperature-dependent,  and  the 
Seebeck  thermoelectric  potential  is  a  function  of  temperature 
gradient.  As  a  result  of  three-dimensional  temperature  distribution, 
the  electric  potential  and  current  density  are  also  three- 
dimensional.  Second,  in  the  energy  equation,  the  source  terms  of 
Joule  heat  and  Thomson  heat  are  the  functions  of  current  density, 
which  are  determined  in  terms  of  the  electric  potential  distribution. 
Thus,  a  multidimensional,  multiphysics  model  with  coupling  of 
temperature  field  and  electric  potential  field  is  needed  to  predict 
the  TE  performance  more  accurately  without  assumption  of 
uniform  current  density  and  to  be  applied  to  various  thermoelectric 
materials  and  comprehensive  operation  conditions. 


Based  on  the  above  analysis,  this  paper  firstly  focuses  on  the 
construction  of  a  general  model  with  coupling  of  temperature  field 
and  electric  potential  field,  which  can  be  used  to  investigate  the 
performances  of  both  the  TEG  and  the  TEC.  Then  a  TEC  is  selected  as 
an  example  to  measure  the  ability  of  the  present  model.  The  effects 
of  the  temperature-dependent  properties  of  p-type  and  n-type 
semiconductors  and  the  convective  boundary  conditions  at  the  TEC 
side  surfaces  on  the  TEC  cooling  capacity  and  coefficient  of 
performance  (COP)  are  analyzed.  The  results  are  compared  to  those 
with  constant  property  assumption  and  adiabatic  boundary 
conditions  extensively  adopted  by  previous  models. 

2.  Model 

A  TEC  device  is  formed  from  a  number  of  thermoelectric 
elements,  connected  electrically  in  series  and  thermally  in  parallel. 
With  consideration  of  periodicity,  a  TEC  element  is  extracted  and  is 
shown  in  Fig.  1,  which  is  composed  of  p-type  and  n-type  semi¬ 
conductor  columns  with  the  same  square  cross-section  sandwiched 
between  two  metallic  interconnectors.  The  interconnectors  and 
semiconductors  have  heights  of  H i  and  H2,  and  the  side  length  of 
the  cross-section  for  semiconductors  is  L2  with  the  distance  of  Lj 
between  p-type  and  n-type  semiconductors. 

The  current  carriers  are  electrons  and  holes  in  the  p-type  and  n- 
type  semiconductors,  respectively.  When  a  direct  current  passes 
from  the  p-type  semiconductor  to  the  n-type  semiconductor,  holes 
in  the  p-type  semiconductors  and  electrons  in  the  n-type  semi¬ 
conductors  are  all  transported  from  bottom  to  top  along  the 
semiconductor  height  direction  (plus  z  axis  direction),  the  corre¬ 
sponding  Peltier  heats  will  be  generated  at  the  interface  between 
interconnectors  and  semiconductors.  The  heat  is  adsorbed  on  the 
bottom  of  semiconductors  and  is  liberated  on  the  top  of  semi¬ 
conductors  due  to  Peltier  effect,  which  causes  that  the  bottom  is 
cooled  and  the  top  is  heated,  and  forms  a  cold  junction  and  a  hot 
junction  with  temperatures  of  TL  and  TH,  respectively.  When  the 
temperature  difference  AT=  Th  -  Ti  is  established  between  the  hot 
and  cold  junctions,  a  Seebeck  electric  potential  is  generated  and  the 
potential  is  proportional  to  the  temperature  difference.  The 
temperature  difference  also  drives  the  heat  to  conduct  back  from 
the  hot  junction  to  the  cold  junction  (Fourier’s  heat  conduction), 
weakening  the  cooling  capacity  of  the  TEC.  The  TEC  element  can  be 
regarded  as  the  combination  of  a  series  of  infinitesimal  Peltier 
junctions  and  each  of  them  individually  absorbs  and  liberates  heat. 
Thus,  under  the  action  of  current  and  temperature  gradient, 
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Fig.  1.  Schematic  diagram  of  the  TEC  element. 
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Thomson  heat  will  be  produced  within  the  TEC  element.  In  addi¬ 
tion,  heat  will  be  generated  along  the  thermoelectric  element  due 
to  Joule  heating,  and  heat  will  also  be  liberated  from  TEC  element  to 
the  ambient  gas  by  thermal  radiation  and  convection.  Thus, 
however,  cooling  capacity  cannot  be  as  strong  as  the  Peltier  heat 
pumping  rate;  it  should  be  determined  by  the  combination  of 
Peltier  heat,  Fourier’s  heat  conduction,  Joule  heat,  Thomson  heat, 
and  heat  losses  to  the  ambient  gas  caused  by  radiation  and 
convection. 

Based  on  the  above  analysis,  a  TEC  model  with  coupling  of 
temperature  field  and  electric  potential  is  constructed  in  the 
present  paper.  It  is  noted  that  the  model  can  also  be  applied  to 
investigate  the  performance  of  the  TEG,  the  difference  from  TEC 
modeling  is  that  heat  should  be  supplied  to  hot  junction  to  produce 
electric  potential  with  the  corresponding  boundary  conditions. 

The  energy  equation  for  the  interconnector,  the  p-type  and  n- 
type  semiconductors  can  be  expressed  as 

I2  -* 

V-(AjVT)  +  —  —  ft J -VT  =  0  (1) 

V-(ApVT)+^-/?p7-VT=  0  (2) 

Op 


for  the  purpose  of  constructing  a  general  model  of  TEC,  the  variable 
properties  should  be  considered.  In  the  present  model,  the  thermal 
conductivity,  electric  conductivity  and  Seebeck  coefficient  are  all 
assumed  to  be  the  functions  of  temperature  as  follows: 

A(T)  =  A(T0)[l+A1(T-T0)+A2(r-T0)2j  (8) 

«  P(T)  =  P(T0)  [l  +  Bi  (T  -  T0)  +  B2  (T  -  ft)2]  (9) 

«(T)  =  a(ft)  [l  +Ct(T - ft)  +  C2(T  -  T0)2]  (10) 

where  To  is  the  reference  temperature;  /.(To),  p(To),  and  a(ft)  are  the 
thermal  conductivity,  electric  resistivity,  and  Seebeck  coefficient  at 
reference  temperature  of  T0,  respectively;  AhA2,  B i,  B2,  ft,  and  C2  are 
the  fitted  parameters  determined  by  experimental  data. 

The  governing  equations  (1)— (3)  and  (5)— (7)  are  solved  itera¬ 
tively  to  obtain  the  distributions  of  temperature  and  electric 
potential  in  the  p-type  and  n-type  semiconductors.  The  boundary 
conditions  adopted  are: 

TiUo  =  Ti  (H) 


V-(AnVT)+^--/?n7-Vr  =  0  (3) 

On 

where  Ai,  Ap,  and  An  are  the  thermal  conductivities  for  inter¬ 
connector,  p-type  and  n-type  semiconductors,  respectively.  The 
second  and  third  terms  on  the  left  side  in  Eqs.  (1)— (3)  denote  Joule 
heat  and  Thomson  heat,  respectively.  J{x,  y,  z)  is  the  local  current 
density,  <7;,  <rp,  and  crn  are  the  electric  conductivities  for  inter¬ 
connector,  p-type  and  n-type  semiconductors,  respectively,  ft,  ft, 
and  ft,  are  the  Thomson  coefficients  for  interconnector,  p-type  and 
n-type  semiconductors  which  can  be  expressed  as: 


(4) 


rilz=2H1+H2  =  ^  (12) 

The  Peltier  heat  will  be  generated  at  the  interface  between 
semiconductor  and  interconnector,  thus,  the  temperature  and  heat 
flux  are  assumed  to  be  continuous  for  this  interface  and  can  be 
expressed  as, 

On  the  interface  of  z  =  Hi,  ft/2  <  x  <  ft/ 2  +  ft,  0  <  y  <  ft: 

1)  \z=Hi  =  TpLh,  (13) 

+ a,Tjziz=H' = ~;'ptei=Hl + “p77ziz=Hi  (i4) 

On  the  interface  of  z  =  Hi,  3ft/2  +  ft  <  x  <  3ft/2  +  2 L2, 

0  <y  <  ft: 


where  a  is  the  Seebeck  coefficient.  It  is  noted  that  constant  Seebeck 
coefficient  implies  that  the  Thomson  effect  is  ignored.  Theoreti¬ 
cally,  the  Thomson  effect  cannot  be  ignored  when  high  tempera¬ 
ture  gradient  occurs  within  semiconductors.  Huang  et  al.  [26]  and 
Chen  et  al.  [43]  have  confirmed  that  TEC  performance  could  be 
improved  not  only  by  increasing  the  figure-of-merit  of  thermo¬ 
electric  materials  but  also  by  taking  the  Thomson  effect  into 
account.  The  transport  of  electrons  and  holes  in  semiconductors  are 
driven  by  electric  potential,  hence,  the  electric  potential  field  is 
needed  to  be  solved  as  follows: 

V-(<r(V0  —  aVT))  =  0  (5) 

where  <p  is  the  electric  potential,  and  aVT  is  Seebeck  electromotive 
force  coming  from  Seebeck  effect.  Once  <j>  is  determined,  the  electric 
field  can  be  calculated  by  the  following  equation: 


z=H,  =  Tn  |z=Hl 


'^Lh,  +“i7Jzlz=H»  =  -AnlFL 

On  the  interface  of  z  =  Hi  +  H2,  ft/2 

Tilz=H,  +H2  =  rp|H,+H2 


r^nl 


+  «n7Jz|z=H, 

;  x  <  Li/2  +  ft,  0  <  y 


(15) 

(16) 


<  ft: 
(17) 


-Aft7 

0zIz=h,+h2 

+  «i77zU„,+H2  (18) 

On  the  interface  of  z  =  Hi  +  H2,  3ft /2  +  ft  <  x  <  3ft/2  +  2 ft, 
0  <  y  <  ft: 


f  “p  #z  Z.  Jf,  I 


Thus,  the  current  density  vector  in  Eqs.  ( 1 )— (3 )  can  be  calculated 
as  follows: 


z=h,+h2 

i.aTnl 


=  /n|, 


+  ®)lJzlz=ff,+H2 


(19) 


(20) 


Minnich  et  al.  [3]  proposed  that  the  properties  of  some  ther¬ 
moelectric  materials  are  strongly  temperature-dependent;  hence 


Volklein  et  al.  [17]  and  Yao  [44]  proposed  that  the  thermal 
radiation  and  convection  heat  transfer  between  the 
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semiconductors  and  the  ambient  gas  cannot  be  negligible  for 
microcooler.  Therefore,  the  convective  boundary  conditions  are 
assumed  on  the  side  surfaces  of  p-type  and  n-type  semiconductors 
and  interconnector  in  the  present  model,  or: 

(21) 

where,  h  is  the  total  heat  transfer  coefficient  including  the  contri¬ 
bution  of  natural  convective  heat  transfer  and  radiation  heat 
transfer,  and  h  =  0  W  m-2  K_1  denotes  adiabatic  boundary  condi¬ 
tions  adopted  in  the  majority  of  previous  models.  Ta  is  the  ambient 
temperature  assumed  to  be  293  K  here. 

The  electric  boundary  conditions  for  TEC  element  are  assumed 
to  be  that  the  inlet  current  is  constant  and  outlet  electric  potential 
is  zero,  or: 

f|x=0  =  constant  0  <y<I2,  0<z<H1  (22) 

</,lx=2i,+2i2  =  0  0  <y  <L2,  0  <  z  <  Hi  (23) 

On  the  other  surfaces,  the  current  cannot  flow  out  of  TEC 
element,  thus,  we  have: 


with  constant  properties  due  to  weak  thermoelectric  effect.  Table  1 
lists  all  the  parameters  in  Eqs.  (8)— (10)  for  the  n-type  and  p-type 
semiconductors  and  the  copper  interconnector. 

3.1.  Grid  independence  examination 

In  order  to  guarantee  that  the  simulation  results  are  indepen¬ 
dent  on  the  grid,  four  different  grid  schemes  were  tested.  The 
geometric  parameters  of  the  TEC  element  and  corresponding  grid 
number  for  test  are  listed  in  Table  2.  The  temperature  difference 
between  the  hot  and  cold  junctions  is  assumed  to  be 
AT  =  Th  -  TL  =  0  K  with  the  total  heat  transfer  coefficient 
h  =  0  W  m-2  K-1.  Fig.  2  shows  the  cooling  capacity  Ql  and  COP  of 
the  TEC  element  as  functions  of  the  current  /  for  four  grid  schemes. 
At  I  <  1.6  A,  four  grid  schemes  have  almost  the  same  Ql  and  COP, 
while  at  I  >  1.6  A,  the  differences  for  Ql  and  COP  predicted  by  four 
grid  schemes  occur.  The  maximum  deviation  between  Grid  (ii)  and 
Grid  (iii)  is  2.18%,  it  reduces  to  0.62%  for  Grid  (iii)  and  Grid  (iv). 
Hence,  to  save  computation  efforts  without  loss  in  accuracy,  the 
Grid  (iii)  is  adopted. 

3.2.  Validation  of  model  self-consistency 


7-t?=0  (24) 

There  are  two  parameters  used  to  evaluate  the  performance  of 
TEC,  cooling  capacity  Ql  and  coefficient  of  performance  (COP).  Ql  is 
defined  as  the  heat  absorbed  from  the  cold  junction,  and  COP  is 
defined  as  the  ratio  of  cooling  capacity  to  electric  power  expressed 


Ql 

W 


(25) 


where  P  is  the  electric  power,  V  is  the  electric  potential  difference 
through  the  TEC  element.  When  side  surfaces  of  TEC  element  are 
adiabatic,  based  on  the  energy  conservation,  we  have: 


The  temperature  filed  and  electric  potential  field  are  coupled  in 
the  present  model,  the  electric  potential  difference  V  through  the 
TEC  element  can  be  solved  directly  as  shown  in  Fig.  3.  Thus,  the 
electric  power  P  can  be  calculated  as  Pi  =  /  x  V.  With  the 
assumption  of  adiabatic  boundary  condition  on  the  side  surfaces  of 
TEC  element,  the  electric  power  can  also  be  calculated  by 
Pj  =  Oh  -  Ql  due  to  energy  balance.  The  electric  powers  calculated 
by  two  approaches  are  compared  and  shown  in  Fig.  3.  The 
maximum  relative  deviation  (P2  -  Pi)/Pi  for  various  currents  is 
only  —0.09%,  indicating  that  the  present  model  has  complete  self- 
consistency. 

3.3.  Comparison  between  numerical  model  and  theoretical  solution 


P  =  Qh-Ql 


(26) 


where  Qh  is  the  heat  dissipated  rate  from  hot  junction.  Thus,  with 
adiabatic  boundary  conditions,  the  COP  is  also  calculated  by  the 
following  equation,  or: 


(27) 


It  is  noted  that  the  COP  was  calculated  only  by  Eq.  (27)  in  the 
previous  models  because  the  electric  potential  difference  is 
unknown. 


In  order  to  verify  the  present  model,  the  temperature  distribu¬ 
tions  predicted  by  the  present  model  and  by  one-dimensional 
theoretical  solution  [26]  are  compared.  In  the  theoretical  solution, 
the  properties  of  semiconductors,  including  the  thermal  conduc¬ 
tivity,  electric  conductivity,  and  Seebeck  coefficient,  are  all  assumed 
to  be  constant  with  adiabatic  boundary  conditions  on  the  side 
surfaces  of  TEC  element.  It  is  noted  again  that  constant  Seebeck 
coefficient  implies  ignored  Thomson  effect,  which  is  applicable  for 
small  temperature  gradient  within  the  semiconductors.  With  these 
conditions,  one-dimensional  theoretical  solution  for  energy  equa¬ 
tion  can  be  expressed  [26]: 


3.  Model  validation 


T(z)-Tl  1  I2H'i  fz  z2\  z 

Th  -  %  *  2  l#2(T h~7l)Vh  HV  H 


(28) 


The  n-type  and  p-type  semiconductors  are  assumed  to  be 
Bi2(Teo.94Seo.o6)3  and  (Bio.25Sbo.75)Te3  with  temperature-dependent 
properties,  respectively.  The  copper  is  selected  as  interconnector 


The  comparison  of  temperature  distributions  between  the 
present  model  and  Eq.  (28)  is  shown  in  Fig.  4,  in  which  the 
geometric  parameters  of  TEC  element  are  the  same  as  Section  3.1 


Table  1 

Parameters  in  Eqs.  (8)-(10)  for  thermoelectric  materials  (T0  =  300  K). 

Material  Thermal  conductivity  Electric  resistivity  Seebeck  coefficient 

_ Aq  (W  ITT1  It"1)  /MK-1)  A2(K-2)  po  (Q  m)  Bt  (K"1)  B2  (1<~2)  tvotVlC1)  C,  (1C1)  C2  (1<-2) 

p-type  1.472  -1.29  x  10-3  1.35  x  10-5  8.826  x  10-6  5.88  x  10-3  8.93  x  10-6  2.207  x  10-4  1.55  x  10-3  -3.15  x  10-6 

n-type  1.643  -9.80  x  10~4  1.56  x  10  5  8.239  x  10  6  4.70  x  10-3  2.67  x  10“6  -2.23  x  10  4  5.62  x  10“4  -4.65  x  10“6 

Interconnector  400  0  0  1.7  x  10~9  0  0  6.5  x  1 0  6  0  0 
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Table  2 

The  geometric  parameters  and  test  grid  numbers  for  TEC. 

Ii  L2  L3  Ht  H2 

Geometry  Value  0.1  mm  0.5  mm  0.2  mm  0.1  mm  1.0  mm 

Grid  i  Grid  number  2  10  4  2  20 

Grid  ii  Grid  number  4  20  8  4  40 

Grid  iii  Grid  number  5  25  10  5  50 

Grid  iv  Grid  number  6  30  12  6  60 


with  Tl  =  300  K,  TH  =  330  K,  h  =  0  W  m  2  K"1,  and  I  =  1.7  A.  The 
temperature  distributions  for  four  different  property  conditions  are 
simulated  by  the  present  model:  1)  all  properties  are  constant  and 
exactly  the  same  as  the  theoretical  solution;  2)  temperature- 
dependent  a  with  constant  A  and  a\  3)  constant  a  with 
temperature-dependent  X  and  rr;  4)  all  properties  are  temperature- 
dependent. 

Fig.  4  shows  that  the  maximum  temperature  occurs  inside 
semiconductors  for  all  four  cases  due  to  thermoelectric  effect.  The 
temperature  of  p-type  semiconductor  is  larger  than  that  of  n-type 
semiconductor  because  Ap  of  p-type  semiconductor  is 
1.472  W  m-1  K-1  at  T  —  300  K,  which  is  lower  than 
An  =  1.643  W  m-1  K-1  for  n-type  semiconductor.  It  can  be  seen  that 
the  temperature  distribution  curves  of  p-type  and  n-type  semi¬ 
conductors  predicted  by  the  present  model  with  the  constant 
properties  and  adiabatic  boundary  conditions  coincide  perfectly 
with  the  theoretical  solution,  validating  the  model.  However,  for 
the  other  three  cases  with  variable  properties,  the  temperatures  are 
higher  than  that  with  constant  properties  and  maximum  temper¬ 
ature  difference  reaches  to  13.25  K.  The  comparison  demonstrates 
that  the  temperature-dependent  properties  should  be  taken  into 
account  to  predict  the  TEC  performance  accurately.  Fig.  4  also 
confirms  the  results  reported  by  previous  investigations  [26,43], 
the  Thomson  effect  could  influence  significantly  on  the  TEC 
performance,  which  decreases  the  temperature  of  semiconductors, 
especially  for  the  p-type  semiconductor  with  lower  thermal 
conductivity. 

4.  Results  and  discussion 

The  cooling  capacity  Ql  and  COP  of  TEC  are  closely  related  to  its 
operating  conditions,  including  the  current  and  the  temperature 
difference  between  the  hot  and  cold  junctions.  The  effect  of  these 
parameters  on  the  TEC  performance  has  been  discussed  extensively 
in  the  previous  investigations,  however,  the  models  do  not  couple 


1(A) 

Fig.  2.  Grid  independence  examinations. 


Fig.  3.  Electric  potential  differences  and  electric  powers  at  various  currents  for  the  TEC 
element. 


the  temperature  field  and  the  electric  potential  field  and  generally 
assumed  constant  properties  of  semiconductors  and  adiabatic 
boundary  conditions  on  the  side  surfaces  of  TEC,  and  hence  the 
results  are  applicable  only  at  specific  conditions.  This  work  will 
analyze  the  effect  of  the  properties  (Section  4.1)  and  convective 
boundary  conditions  (Section  4.2)  on  the  performance  of  the  TEC 


Fig.  4.  Comparison  of  temperature  distributions  predicted  by  the  present  model  and 
calculated  by  the  one-dimensional  theoretical  solution:  (a)  p-type  semiconductor;  (b) 
n-type  semiconductor. 
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element  in  a  wide  range  of  the  current  and  temperature  difference 
between  the  hot  and  cold  junctions  using  the  coupled  model 
described  in  Section  2. 

4.7.  Effect  of  semiconductor  properties 

The  three-dimensional  analysis  is  carried  out  with  TEC  element 
listed  in  Table  2.  The  temperature  of  cold  junction  is  fixed  to 
Tl  =  300  K  and  the  temperature  difference  between  the  hot  and 
cold  junctions  is  assumed  to  be  AT  =  0  K,  20  K,  40  K,  60  K,  80  K,  and 
100  K,  respectively.  The  adiabatic  boundary  condition, 
h  =  0  W  m-2  K_1,  is  adopted  in  this  section  to  highlight  the  effect  of 
semiconductor  properties. 

Fig.  5  illustrates  the  cooling  capacities  of  the  TEC  element  as 
a  function  of  current  with  various  AT  under  constant  and  variable 
property  models,  respectively.  For  the  constant  property  model,  Ql 
firstly  increases  and  then  drops  as  the  current  further  increases, 
hence  there  exists  maximum  cooling  capacity  Q^max-  At  AT  -  0  K, 
effective  working  currents,  defined  as  the  current  when  Ql  >  0  W, 
are  0  A  <  /  <  3.940  A,  and  the  current  corresponding  to  QL,max  is 
(3.940  -  0)/2  =  1.970  A,  implying  that  the  Ql— 7  curve  is  completely 
symmetric.  As  AT  is  increased,  the  effective  working  current 
interval  is  reduced  but  the  Ql—I  curve  remains  symmetric  with  the 
same  current  of  1.970  A  corresponding  to  Qi^max-  The  maximum 
cooling  capacity  QL,max  is  reduced  as  AT  is  increased  due  to  elevated 
Fourier’s  heat  conduction.  The  constant  decrease  rate  (dQL,max/dT) 
of  7.75  x  10-4  W/K  is  observed,  and  similar  results  were  also  re¬ 
ported  by  Lee  and  Kim  [30],  For  example,  at  AT  =  0  K, 
0L,max  =  0.1302  W;  at  AT  =  20  K, 

0L,rnax  =  0.1302  -  7.75  x  10“4  x  20  =  0.1147  W;  at  AT  =  100  K, 
Qcmax  =  0.1302  -  7.75  x  10“4  x  100  =  0.0527  W. 

Fig.  5  indicates  that  the  Ql—I  curves  for  variable  property  model 
are  distinct  from  those  for  constant  property  model.  First,  the 
maximum  effective  working  current  is  decreased  significantly,  the 
Ql— 7  curves  no  longer  have  symmetry.  At  low  currents,  the  Joule 
heat  and  Thomson  heat  produced  within  the  semiconductors  are  so 
small  that  the  temperature  distributions  predicted  by  the  constant 
and  variable  property  models  are  nearly  identical  (Fig.  6),  leading  to 
a  small  difference  in  cooling  capacity  Ql  between  the  two  models. 
As  the  current  is  raised,  the  temperature  distributions  predicted  by 
the  two  models  become  remarkably  different  because  the  Joule 
heat  and  Thomson  heat  are  increased  rapidly,  and  the  variable 
property  model  predicts  a  higher  temperature  distribution.  Higher 
temperature  distribution  results  in  a  stronger  Fourier’s  heat 


Fig.  5.  Cooling  capacity  of  the  TEC  element  for  constant  and  variable  properties  at 
various  temperature  differences. 


conduction  from  hot  junction  to  cold  junction,  thus,  the  cooling 
capacity  Ql  and  the  corresponding  maximum  effective  working 
current  are  decreased.  Second,  maximum  cooling  capacity  Qi^max  is 
reduced  significantly  for  the  variable  property  model,  and  the 
currents  corresponding  to  Q^max  at  different  AT  are  no  longer  the 


a 


temperature  differences:  (a)  constant  properties;  (b)  variable  properties. 
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same  but  decrease  as  AT  increased.  For  example,  at  AT  =  0  K, 
QL,max  =  0.1220  W,  the  corresponding  current  is  1.685  A;  at 
AT  =  100  K,  QL,max  —  0.0431  W,  the  corresponding  current  decreases 
to  1.575  A.  Furthermore,  similar  as  the  constant  property  model, 
0L,max  is  decreased  as  AT  increased  and  the  decrease  rate  still 
remains  constant  value  of  7.89  x  10-4  W/I<  slightly  lower  than 
7.75  x  10”4  W/K  for  the  constant  property  model. 

Fig.  7  shows  the  COPs  of  TEC  element  as  functions  of  current  and 
temperature  difference  for  the  two  models.  At  AT  =  0  K,  the  COP  is 
reduced  monotonously  as  the  current  increased  for  both  models, 
and  there  is  no  maximum  COP.  However,  at  AT  >  0  K,  the  COP 
reaches  its  peak  value  COPmax  at  a  certain  current  and  then 
monotonically  drops  as  the  current  increases  further.  It  is  also  noted 
that  the  COP  decreases  as  AT  increased.  The  I—V  curves  of  TEC 
element  for  the  two  models  are  shown  in  Fig.  8.  It  is  seen  that  the 
electric  potential  difference  through  the  TEC  element  increases  as 
AT  increased  for  two  models.  Thus,  higher  AT  not  only  decreases 
the  cooling  capacity  but  also  increases  the  electric  potential 
difference  leading  to  the  reduced  COP.  Fig.  8  also  shows  that  the 
electric  potential  difference  increased  almost  linearly  with  the 
current  for  constant  property  model,  which  is  lower  than  that  for 
variable  property  model,  the  discrepancy  between  two  models 
becomes  more  significant  at  high  currents  and  large  AT.  The 
discrepancy  of  the  electric  potential  differences  comes  from  the 
difference  of  temperature  distributions  predicted  by  two  models. 
The  Seebeck  electric  potential  gradient  is  equal  to  the  product  of 
Seebeck  coefficient  and  the  temperature  gradient.  The  higher 
temperature  gradient  and  positive  temperature  effect  of  Seebeck 
coefficient  cause  a  larger  Seebeck  potential  through  the  TEC 
element  for  the  variable  property  model.  It  is  noted  from  Fig.  7  that 
the  current  corresponding  to  COPmax  is  elevated  as  AT  increased, 


and  this  current  is  different  from  the  current  corresponding  to 
0L,max.  which  means  that  the  requirement  of  maximum  cooling 
capacity  and  COP  cannot  be  met  simultaneously. 

It  can  be  seen  from  Fig.  7  that  for  AT  =  20  K,  COPmax  =  2.403  at 
current  of  0.270  A  for  the  variable  property  model,  and 
COPmax  =  2.434  at  current  of  0.281  A  for  the  constant  property 
model  with  relative  deviation  of  (2.434  -  2.403)/2.403  =  1.3%. 
However,  at  AT  =  100  K,  COPmax  =  0.187  at  current  of  1.116  A  for  the 
variable  property  model,  and  COPmax  =  0.220  at  current  of  1.262  A 
for  the  constant  property  model  with  relative  deviation  of 
(0.220  -  0.187)/0.187  =  17.6%.  Therefore,  the  constant  property 
model  overestimates  the  COPmax  and  the  degree  of  the  over¬ 
estimation  increases  with  AT  increased.  This  can  be  explained  by 
the  following  two  reasons.  First,  as  AT  increases,  the  constant 
property  model  overestimates  the  cooling  capacity  (Fig.  5)  and 
underestimates  the  electric  potential  difference  (Fig.  8)  more 
significantly.  Second,  the  current  corresponding  to  COPmax 
increases  at  larger  AT,  higher  current  further  increases  the  degree  of 
overestimation  of  the  constant  property  model.  The  above  results 
indicate  that  the  temperature-dependent  properties  of  semi¬ 
conductors  are  necessary  to  take  into  account  for  TEC  modeling  at 
large  AT  and  higher  I. 

4.2.  Effect  of  convective  boundary  condition 

During  TEC  operation,  the  temperature  of  TEC  is  generally 
higher  than  the  ambient  temperature;  hence,  the  heat  will  be 
transferred  from  the  TEC  to  the  ambient  gas  by  radiation  and 
convection.  In  the  previous  models,  the  adiabatic  boundary 
conditions  were  adopted  extensively,  only  a  few  models  considered 
the  effect  of  radiation  and  convection  heat  transfer  on  the  TEC 
performance  [17,26,38],  In  these  models,  the  radiation  heat  transfer 
hr  was  linearized  under  the  assumption  of  small  temperature 
difference  ((TPin  -  Too  )/T „  «  1 )  and  a  total  heat  transfer  coefficient 
combining  the  radiation  heat  transfer  coefficient  hr  and  the 
convection  heat  transfer  coefficient  hc,  which  can  be  expressed  as 
h  =  hc  +  hr  =  hc  +  4e<7BT„,  is  used  to  estimate  the  heat  transfer 
between  the  TEC  and  the  ambient.  However,  different  values  of  h 
were  assumed  in  these  models,  for  example,  h  is  assumed  to 
belOO  W  m~2  K_1  in  Ref.  [17],  50  W  m"2  K"1  in  Ref.  [26],  and  5- 
15  W  m-2  K_1  in  Ref.  [38].  In  the  present  model,  the  value  of  h  is 
assumed  to  be  varied  from  0  to  200  W  m-2  K-1,  it  is  noted  that 
a  little  higher  h  (=200  W  m-2  K-1)  is  employed  here  to  highlight 
the  effect  of  convective  boundary  condition. 

The  calculations  are  carried  out  with  the  geometry  of  TEC 
element  identical  to  that  in  Section  4.1.  The  temperature  of  the  cold 
junction  is  set  at  300  K,  and  the  temperature  of  the  hold  junction  is 
300  K,  340  K,  and  400  K,  respectively.  Thus,  the  temperature 
difference  is  0  K,  40  K,  and  100  K,  respectively.  In  addition,  the 
temperature-dependent  properties  are  adopted. 


Fig.  9.  Effect  of  the  total  heat  transfer  coefficient  on  the  temperature  distributions  on  the  middle  cross-section  (y  =  0.25  mm)  at  AT  =  0  K  and  I  =  0.8  A:  (a)  h  =  0  W  m-2  K-1;  (b) 
h  =  50  W  m-2  IT1;  (c)  h  =  200  W  m~2  K“l. 
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a  b  c 


Fig.  10.  Effect  of  the  total  heat  transfer  coefficient  on  the  temperature  distributions  on  the  middle  cross-section  (y  =  0.25  mm)  at  AT  =  0  K  and  1  =  2.2  A:  (a)  h  =  0  W  m  2  K  (b) 
h  =  50  W  nr2  K  (c)  h  =  200  W  mr2  K-1. 


Figs.  9  and  10  show  the  effect  of  the  total  heat  transfer  coeffi¬ 
cient  on  the  temperature  distributions  on  the  middle  cross-section 
( y  =  0.25  mm)  at  AT- OK  with  1  =  0.8  A  and  2.2  A,  respectively.  It  is 
seen  again  that  the  p-type  semiconductor  has  higher  temperature 
than  the  n-type  semiconductor  due  to  its  lower  thermal  conduc¬ 
tivity  and  higher  electric  resistivity.  At  low  current  of  /  =  0.8  A,  there 
is  a  small  amount  of  Joule  heat  and  Thomson  heat  produced  inside 
semiconductors,  the  difference  of  temperature  distributions  for 
various  total  heat  transfer  coefficients  is  small.  As  the  total  heat 
transfer  coefficient  increases,  the  temperatures  of  the  semi¬ 
conductors  are  reduced  slightly  since  the  heat  loss  to  the  ambient  is 
enhanced.  With  h  =  0  W  m-2  K_1,  the  maximum  temperature  for 
the  p-type  semiconductor  is  Tpmax  =  308.0  K  and  for  the  n-type 
semiconductor  is  Tnmax  =  306.6  K;  with  h  =  50  W  m-2  K-1, 
Tp,max  =  307.6  K,  and  Tn,max  =  306.3  K;  with  h  -  200  W  m"2  K-1, 
Tp’max  =  306.4  K,  Tn,max  =  305.3  K.  At  high  current  of  I  =  2.2  A, 
however,  the  Joule  heat  and  Thomson  heat  produced  inside  the 
semiconductors  increase  significantly,  causing  that  the  tempera¬ 
tures  of  both  the  p-type  and  n-type  semiconductors  are  far  larger 
than  those  at  I  =  0.8  A.  Increased  temperature  difference  between 
the  semiconductors  and  the  ambient  enhances  the  heat  loss  to  the 
ambient  and  hence  increases  the  effect  of  total  heat  transfer  coef¬ 
ficient.  The  above  analysis  can  be  proved  by  the  fact  that  with 
h  =  0  W  m”2  K"1,  Tpjnax  =  385.8  K,  and  Tn,max  =  361.2  K;  with 
h  =  50  W  m  2  K-t  Tp,max  =  382.3  K,  and  Tn,max  =  359.2  K; 
h  =  200  W  m  2  K"1,  Tp,max  =  373.1  K,  and  Tn,max  =  354.0  K.  In 
addition,  it  is  observed  that  at  high  total  heat  transfer  coefficient  of 
h  =  200  W  m-2  K-1,  the  isothermal  lines  near  the  middle  of 


semiconductors  are  not  parallel  to  the  x-direction,  indicating 
a  three-dimensional  temperature  distribution  occurs. 

Fig.  11  shows  the  I—V  curves  of  the  TEC  element  for  various  h  at 
AT  =  0  K.  At  /  <  1.5  A,  h  has  a  very  small  effect  on  the  I—V  curves, 
while  at  /  >  1.5  A,  the  effect  of  h  is  elevated  gradually  as  the  current 
increases.  Higher  h  leads  to  a  reduced  electric  potential  difference 
through  the  TEC  element  and  hence  a  reduced  electric  power. 

The  cooling  capacity  Ql  and  the  corresponding  COP  for  various  h 
at  AT=  0 1<  are  shown  in  Fig.  12.  Similar  as  the  electric  potential,  the 
effect  of  h  on  Ql  is  enhanced  with  the  current  increased.  Increase  in 
h  leads  to  higher  Ql  because  the  heat  loss  to  the  ambient  decreases 
the  Fourier’s  heat  conduction  to  cold  junction.  Fig.  12  also  indicates 
that  as  h  is  raised,  the  effective  working  current  interval  increases 
and  the  current  corresponding  to  QL,max  moves  towards  to  higher 
current  direction.  At  h  =  0  W  m-2  K-1,  the  interval  is  varied  from 
0  A  to  2.806  A  with  QL,max  =  0.1220  W  at  the  current  of  1.685  A;  at 
h=  50  W  nr2  K  the  interval  is  0-2.858  A  with  Oi  =  0.1256  W 
at  the  current  of  1.722  A;  at  h  =  200  W  m-2  I<-1,  the  interval  is  0— 
3.014  A  with  QL,max  —  0.1355  W  at  the  current  of  1823  A.  Fig.  12 
demonstrates  that  the  COP  increases  with  h  increased  due  to  that 
higher  h  not  only  increases  Ql  but  also  decreases  V. 

The  effect  of  h  on  the  TEC  performance  is  closely  dependent  on 
AT.  Figs.  13  and  14  show  the  Ql  and  the  corresponding  COP  for 
various  h  at  AT  =  40  K  and  AT  =  100  K,  respectively.  Similar  as 
AT  =  0  K,  Ql  and  COP  both  increase  as  h  increases  and  the  effect  of  h 
becomes  more  significant  at  higher  currents.  The  above  analyses 
indicate  that  the  radiation  and  convection  heat  transfer  between 
the  TEC  and  the  ambient  gas  are  beneficial  to  improve  the  TEC 
performance.  It  can  be  expected  that  the  effect  of  the  radiation  and 
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Fig.  13.  Performance  of  the  TEC  element  with  various  convective  heat  coefficients  at 
AT  =  40  K. 


Fig.  14.  Performance  of  the  TEC  element  with  various  convective  heat  coefficients  at 
AT  =  100  K. 


element  with  the  ratio  of  cross-sectional  area  to  the  height  of 
0.25  mm  is  investigated  and  Bi2(Teo.94Seo.o6)3  and  (Bio.25Sbo.7s)Te3 
are  selected  as  the  n-type  and  p-type  thermoelectric  materials, 
respectively. 

The  temperature  filed  and  electric  potential  field  are  calculated 
coupled  in  the  present  model,  so  the  electric  potential  difference 
through  the  TEC  element  can  be  solved  directly.  The  self- 
consistency  of  the  model  was  verified  by  the  fact  that  the  electric 
powers  calculated  by  two  approaches  are  identical  under  adiabatic 
boundary  conditions.  With  constant  properties  and  adiabatic 
boundary  conditions,  the  temperature  distributions  within  the  p- 
type  and  n-type  semiconductors  predicted  by  the  present  model 
agree  well  with  the  one-dimensional  theoretical  solution. 

The  effects  of  the  properties  and  convective  boundary  condi¬ 
tions  on  the  TEC  performance  are  discussed  in  a  wide  current  and 
temperature  difference  range  using  the  present  model.  The  pre¬ 
dicted  results  show  that  the  constant  property  model  over¬ 
estimates  the  cooling  capacity,  the  COP,  and  the  effective  working 
current  interval  of  the  TEC  and  underestimates  the  electric  poten¬ 
tial  compared  to  the  variable  property  model.  Besides,  a  signifi¬ 
cantly  increased  deviation  is  observed  at  high  currents.  The 
convective  boundary  conditions  cause  a  three-dimensional 
temperature  distribution  within  semiconductors.  The  effect  of  the 
radiation  and  convection  heat  transfer  on  the  TEC  performance 
increases  significantly  at  high  temperature  difference  between  the 
hot  and  cold  junctions  or  at  high  current. 
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convection  heat  transfer  for  the  TEG  is  opposite  to  the  TEC,  because 
the  heat  loss  to  the  ambient  gas  reduces  the  temperature  difference 
between  the  hot  and  cold  junctions  and  hence  reduces  the  Seebeck 
electric  potential.  The  comparison  of  Figs.  12—14  indicates  that  the 
effect  of  h  on  the  TEC  performance  is  elevated  as  AT  increases.  The 
reason  can  be  explained  as  follows.  As  mentioned  before,  at  high  AT, 
there  is  a  stronger  Fourier’s  heat  conduction,  which  causes  the  heat 
transferred  back  from  the  hot  junction  to  the  cold  junction  and 
hence  leads  to  a  decreased  Ql.  With  the  same  h,  more  heat  is 
absorbed  by  the  ambient  gas  at  higher  AT  due  to  lager  temperature 
difference  between  the  TEC  and  the  ambient  gas,  so  the  effect  of 
Fourier’s  heat  conduction  is  reduced  and  Ql  is  increased  more 
significantly.  Figs.  13  and  14  again  confirm  that  increase  in  h 
increases  the  effective  working  current  interval  and  causes  a  move 
of  the  current  corresponding  to  Qi^max  towards  higher  current. 
Oppositely,  increase  in  h  causes  that  the  current  corresponding 
COPmax  moves  towards  lower  current. 


5.  Conclusions 

This  paper  develops  a  general,  three-dimensional  TE  model  with 
coupling  of  the  temperature  field  and  the  electric  potential  field.  All 
effect  occurred  in  the  TE  are  taken  into  account,  such  as  the  Joule 
heating,  the  Thomson  effect,  Fourier’s  heat  conduction,  the  Peltier 
effect,  the  radiation  and  convection  heat  transfer,  and  temperature- 
dependent  properties  of  semiconductor  materials.  A  miniature  TEC 
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